Dynamics of doublon-holon pairs in Hubbard two-leg ladders 
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The dynamics of holon-doublon pairs is studied in Hubbard two-leg ladders using the time- 
dependent Density Matrix Renormalization Group method. We find that the geometry of the 
two-leg ladder, that is qualitatively different from a one-dimensional chain due to the presence of a 
spin-gap, strongly affects the propagation of a doublon-holon pair. Two distinct regimes are identi- 
fied. For weak inter-leg coupling, the results are qualitatively similar to the case of the propagation 
previously reported in Hubbard chains, with only a renormalization of parameters. More interesting 
is the case of strong inter-leg coupling where substantial differences arise, particularly regarding the 
double occupancy and properties of the excitations such as the doublon speed. Our results suggest 
a connection between the presence of a spin gap and qualitative changes in the doublon speed, 
indicating a weak coupling between the doublon to magnetic excitations. 

PACS numbers: 71.10.Fd, 71.10.Li, 71.35.-y 



I. INTRODUCTION 

Quasi-one dimensional correlated quantum systems 
have received considerable attention from the Condensed 
Matter community because they display exotic properties 
that arise from the dimensional confinement of electrons. 
A famous example are the quantum spin systems with 
the geometry of "ladders", namely infinitely long one- 
dimensional "legs" that are coupled along "rungs" via 
tight binding hopping terms with strengths comparable 
or larger than those along the legs Previous research has 
shown that ladders with an odd number of legs behave 
similarly as truly one-dimensional systems in the sense 
that spin-spin correlations decay with distance following 
power laws. However, for even number of legs (two in 
particular) the decay i s ex ponential due to the presence 
of a so-called spin gapP^l The existence of this spin gap 
has been confirmed experimentally in materials such as 
the copper oxide SrCu203. For the case of two-leg lad- 
ders, the spin gap is caused by the dominance of spin 
singlet configurations along the rungs, which have a spin 
gap between singlet and triplet states. The existence of 
a ground state dominated by spin singlet arrangements 
resembling resonant valence bond states, and other prop- 
erties such as the prediction of superconductivity upon 
doping, have established interesting analogies between 
two-leg ladders and the two-dimensional high temper- 
ature superconductors based on Cu02 layersP Two-leg 
ladders can also be found in the context of organic com- 
pounds such as BPCB ((C 5 U 12 N) 2 Cu Br 4 )P showing 
that the importance of two-leg ladders is not only re- 
stricted by its similarity with the high-T c cuprates but it 
is broader, covering several families of materials. 

The studies of two-leg ladders are not restricted to spin 
systems, but considerable work has been devoted to the 
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FIG. 1: (Color online) Creation of a doublon-holon pair in a 
two-leg Hubbard ladder. 



exploration of both the spin and the charge degrees of 
freedom in the context of Hubbard models defined on 
two-leg ladder geometries. In fact, computational results 
for Hubbard ladders with strong rung (or inter-leg) in- 
teraction show the presence of a spin gap for all inter- 
action strengthsP^ However, not much is known about 
the real time dynamics of charge excitations in ladder 
systems. Given the prospect that highly correlated ox- 
ides can become technologically useful for light-to-energy 
conversion, the understanding of this interplay is very im- 
portant. A related issue is the role of charge excitations 
in the transport properties of strongly correlated nanos- 
tructures, as highlighted by recent experiments!^ 

A crucial question is whether charge excitations in 
a correlated insulator will be able to properly transfer 
the charge into the metallic contacts, thus establishing 
a steady-state photocurrent. This has been a motivat- 
ing factor for studies of exciton-like pairs of charge ex- 
citations in Mott insulators E2HH1 j n Hubbard systems, 
the elementary charge excitation that preserves the total 
number of electrons is a "doublon-holon" pair, formed 
by a double-occupied site (doublon) and an "empty" one 
(holon) , as shown schematically in Fig. [I] In some sit- 
uations, the holon and doublon can attract each other 
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forming a "Hubbard exciton" , as revealed through spec- 
troscopic evidence by a recent experiment in the transi- 
tion metal oxide YV0 3 P 

Thus, it is important to understand how interactions, 
quantum fluctuations, and dimensionality can affect the 
lifetime of such charge excitations. In our previous 
effort ESI it was shown that the mechanism for holon- 
doublon decay into magnetic excitations is inefficient in 
strictly one-dimensional Hubbard chains. In the present 
paper, we extend these ideas to investigate the dynamics 
of doublon-holon pairs in Hubbard two-leg ladders, which 
are qualitatively different from chains due to the presence 
of a spin gap. We find that the geometry of the ladder 
strongly affects the propagation of a doublon-holon pair 
and identify two distinct parameter regimes. 

For weak inter-leg coupling, the re sults are qualita- 
tively similar to the case of the chairP^lUJ with some 
renormalization of the parameters. More interestingly 
is the case of strong inter-leg coupling: in this regime, 
that is drastically different from a one-dimensional chain, 
quantities such as the double occupancy and the proper- 
ties of the excitations doublon (e.g., the doublon speed) 
are qualitatively affected by the ladder geometry. 

The paper is organized as follows: the Hubbard ladder 
model and the specifics technical details of the tDMRG 
method used are presented in Sec. |TT] and in Appendix 
|A") In Sec. |III[ we present the results for the charge 



(Sec. Ill A) and double occupancy dynamics (Sec. IIIB) 
as well as the dependency of the excitation speed with 
the spin-gap model (Sec. IIIC). Our concluding remarks 
are presented in Sec. |IV) 



II. MODEL AND METHODS 

The Hubbard ladder Hamiltonian is written as 

where the notation is similar to that of Ref. [TUJ The 
hopping matrix t, however, now corresponds to that of a 
two-leg ladder, such that tij = t x if i and j are nearest 
neighbors along the x direction, t^ = t y if i and j are 
nearest neighbors along the y direction, and (zero) oth- 
erwise. In the following, we set t x = 1 and use it as our 
energy unit. 

In order to create charge excitations, we define holon- 
and doublon-creation operators as h\ — ^ a 6^(1 — hia), 
and d\ = J2 a ^ia^is respectively, where f As in 

Ref. HOI we model the electron- and hole-like excitations 
(created by light absorption, for instance) as an excited 
state l^e) = ?ijd]|^o)) where l^o) is the ground state 
of the Hubbard-ladder system at half-filling, and i and 
j are chosen fixed sites where the excitation occurs. For 
concreteness, as depicted in Figure [l] we choose these 
sites to be the central sites of the upper leg, which allows 
a direct comparison to the results depicted in Refs. llOlTTI 
for the case t y = 0. 



Next comes the time evolution of the excited state with 
the tDMRG method!- 14 * 15 ^ In order to time-evolve the sys- 
tem, we will use the time-step-targetting Krylov method 
(Ref. [T5] and references therein) and follow the imple- 
mentation found in Ref. 1171 We find this approach to 
be better suitable for generic Hamiltonians in the ladder 
geometry rather than the Suzuki- Trotter decomposition 
method, which would require the ladder to be treated as 
a series of coupled dimersP We have also compared some 
results to a Runge-Kutta approach and found excellent 
agreement. 

We calculate the state |*(r)) = e" lffT |* e ), with H 
given by Eq. (n|, as a function of time r (such that 
|\I/(0 + )) = j^e)) and then compute real observables as 
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shown are for a 20-site (10 x 2) ladder at half-filling, 
with time step At — 0.2 and times up to r < 6 (in units 
of h/t x ). As we discuss in Appendix |Aj accumulated 
errors during time evolution introduce a constraint in 
the accuracy of the results at long time scales. For this 
reason, we have opted to show results with times only up 
to r max = 6 (see Appendix [A| for details). 



III. RESULTS 

In this section, we show results for the local charge and 
double occupancy at each site, defined respectively as 

(niKr) = J2 (*(T)\cLc ta \y(T))/(*(T)\*(T)) , 

<r=t4 

(nt)(r) = (*(r)|n it n i4 |*(r)>/<*(r)|*(r)) , (2) 

using the state |*(r)) = e~ iHT \^(0 + )). The calculations 
are performed for a 20-site (10 x 2) ladder at half-filling 
unless otherwise noted. In our convention, the sites are 
labeled as follows: the leftmost upper site is site 0, the 
leftmost lower site is site 1, and the site index increases 
by one in the vertical top-down direction and by two in 
the horizontal-right direction. Thus, even numbered sites 
(starting at "0") are located at the upper leg and odd- 
numbered sites are in the lower leg. The doublon-holon 
pair was created at time r = + in the center sites of 
the upper leg, indicated as sites "8" (holon) and "10" 
(doublon). 



A. Charge dynamics 

The typical behavior of the charge (nj)(r) as a function 
of time is shown in Fig. [2] for U — 10. Each panel shows 
the charge configuration of the ladder at a given time 
r for t y = (left panels) and t y = 3 (right panels). 
The top panels show the creation of the doublon ((n) = 
2) and holon ((n) = 0) excitations in the central sites 
of the upper leg at r = + . As time progresses, the 
excitations propagate throughout the system in opposite 



3 



1.5 



(a) 






[(d) 














\1 










4 8 12 16 
Top site i 



4 8 12 16 
Top site i 



FIG. 2: (Color online) Charge (n 4 ) at each site for different 
times for the cases t y = (a,b,c) and t y = 3 (d,e,f). Charge 
values in the upper leg are shifted by 0.2 for better visualiza- 
tion. Times are r = (panels a,d), r = 1 (b,e), and r = 2 
(c,f). 



directions,^ eventually reaching the end of the ladder at 
r~3-4. 

More instructive is to show the average charge as a 
function of time r on a particular site. Fig. [3] shows 
(nj(r) at the sites close to the holon-doublon creation 
point for U = 10, t x = 1 and t y = (corresponding to 
the case of uncoupled ladder legs, with zero spin-gap) 
and t y — 3 (spin gapped case). 

For t y = (left panels in Fig. [3| and r = 0, the 
charge at holon and doublon sites is and 2, respec- 
tively with (rij) = 1 in the remaining sites. As time 
evolves, the charge on the doublon (holon) site decreases 
(increases) and oscillates as time progresses. More impor- 
tantly, the excess (missing) charge is initially transferred 
to the neighboring sites in the same leg (labeled "6" and 
"12") so that the charge excitations "move" throughout 
the upper leg and eventually, the local charge on all sites 
in the upper leg will fluctuate, while keeping the total 
charge N = Yj%{ n i) constant. As expected, for t y = 
there is no transfer to the lower leg (Fig[3~l 



For t y ^ 0, the charge excitations will also propagate 
into the lower leg (Fig. |3lc,d). In the case t y = 3t x , the 
dynamics of the local charge on the sites neighboring, 
say, the doublon site (site "10") vertically or horizontally 
will be different, with a large portion of charge oscillating 
vertically due to the dominance of the hoppings along the 
y direction (t y > t x ). This can be seen in Fig. |3jc,d) by 
comparing the charge on sites "11" (first neighbor in the 
lower leg) and "12" (first neighbor in the upper leg). 

We should also point out that, since the time-evolution 
conserves energy, there is no mechanism for "recombina- 
tion" of the doublon- holon pair once it is created.^ Thus, 
the electron- and hole-like excitations do not "cross" each 
other, creating a spatial particle-hole asymmetry in the 
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FIG. 3: (Color online) Charge (m) on sites close to the 
holon-doublon creation sites (8 and 10) for t y — (a,b), and 
ty = 3 (c,d). 



system. This asymmetry is kept throughout the time evo- 
lution such that the charge of sites in the "electron side" 
are related to their "hole side" counterparts by a particle- 
hole transformation (e.g., (n i= i 2 )(T) = 2 — (n i=6 )("r)). 



B. Double Occupancy 

The behavior of the double occupancy (nf) = (n^n^) 
is shown in Figs. |4]and[5| At r = + , the double occu- 
pancy is the maximum (= 1) at the doublon site, zero 
at the holon site and small in the rest of the system (the 
actual value will depend on U). For t y = 0, the local 
double occupancy fluctuates as the doublon-holon pair 
propagates in the upper leg but there is no change in the 
double occupancy in the lower leg. This is in contrast 
with the behavior for t y > t x , for which the propaga- 
tion of the doublon in the upper leg changes the double 
occupancy of the sites in the lower leg. 

An interesting property seen for t y ^ is that the 
vertical coupling correlates the double occupancy of sites 
within the same "rung" . This is clearly shown in Fig.[4]-d, 
which shows the double occupancy at r = + (compare it 
with the t y = case Fig. |4]-a): the creation of a holon or 
a doublon on a given site also changes the double occu- 
pancy in the site directly below it. This happens even in 
the non-interacting case and it is a natural consequence 
of the dominance of the rungs in the t y 3> t x regime. 

This effect can be better visualized by looking at the 
double occupancy for individual sites as a function of 
time (Fig. [fj). For t y = 3 (Fig. [5}c,d) and time r = 0+, 
when the holon-doublon pair is created at sites 8 and 10, 
the double occupancy at sites in the same rungs (9 and 
11, respectively) is reduced to zero. This is in contrast 
with the other sites in the lower leg (say, sites 7 and 13), 
which retain the ground-state double occupancy value. 
This is quite different from the case of uncoupled legs 



4 




12 
Top site i 



4 8 12 
Top site i 




2 3 

time x 



FIG. 5: (Color online) Double occupancy versus time at sites 



FIG. 4: (Color online) Double occupancy (n^n^) at each close to the holon-doublon creation sites (8 and 10) for t y 







site for different times for the cases t y = (a,b,c) and t y = 3 
(d,e,f). Double occupancy values in the upper leg are shifted 
by 0.2 for better visualization. Times are r = 0(a,d), l(b,e), 
and 2(c,f). 



(t y = 0, Fig. [5]-a,b), for which no change in the double 
occupancy happens in the lower leg, as expected. 

Another interesting feature of the large t y regime is the 
"beating" of both charge and double occupancy within 
the rung where the doublon is created. Figs. [3]-(c,d) and 
[5]-(c,d) show that, as time evolves, the charge and double 
occupancy in site 11 oscillate with the same period but 
out of phase with their counterparts in site 10. Interest- 
ingly, for the "holon rung" (sites 8 and 9) these beatings 
show up only in the charge (Fig. [3}c,d) but not in the 
double occupancy (Fig. |5[c,d). 

Next, we look at the total double occupancy A^(r) = 
J2i( n t)( T )- I n contrast with the total charge (which is 
conserved and therefore independent of r), Nd(r) fluc- 
tuates with time r for U ^ 0. As in Ref. QUI we 
look at the quantity AN d (r) = N d (r) - Nf ] where 

N^' = ($o\n>f l^o) is the total ground-state double oc- 
cupancy. This "excess double occupancy" is shown in 
Fig. [6] 

We note that the value of ANd(0) will depend on 
t y . For t y — 0, the creation of the holon-doublon pair 
only changes the double occupancy locally and therefore 
AN d (0) = l-((*oKJ*o) + (*o|n? d |*o» where ih and 
id label the sites where the holon and doublon are cre- 
ated, respectively. By contrast, for t y ^ 0, the creation of 
doublon-holon pair will change the double occupancy on 
the other sites of the rungs as well (as discussed above) 
by an amount that increases with t v . Thus, ANd(0) de- 
creases with t y , as shown in the inset in Fig. |6| The 
decrease is much accentuated for t y > 2, establishing 
a clear qualitative distinction from the results for the 
one-dimensional case t y = 0, where ANd(0) is "local" in 
the sense that it depends only on the ground-state dou- 
ble occupancy at the sites where the holon-doublon is 
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created.^ 

For longer times, 
height that also decreases with t y . For U — 10 and large 
t y , A-ZV^t) fluctuates around a value that is essentially 
independent of r, with no sizeable decay in double oc- 
cupancy with time. This is shown in the inset of Fig. [6] 
where ANd(r = + ) and AiV £ j(r = 6) are plotted as a 
function of t y . For t y = and different values of U, the 
value of (ANd) T increases with U, consistently with the 
findings of Ref. [lOj 



C. Doublon speed 



A clear difference in behavior between the two regimes 
(small and large t y ) can be captured by the speed at 
which the charge excitations travel, particularly the dou- 
blon. This is illustrated in Fig. [7] which shows the double 
occupancy vs time at site "18" , located in the upper leg 
and at the right end of the ladder. 

The position of the first peak in {na)(r) marks the "ar- 
rival time" T a of the doublon and can be used to infer its 
"speed" Vd car" 1 . In Fig.j7[ we can identify two distinct 
regimes: for t y < 1 (Fig. [Wa), such arrival times are of 
order r a ~ 3 while for t y > 1 (Fig. [7]-b), arrival times are 
substantially higher, in the 3 < r a < 4 range. Thus, the 
doublon effectively moves "faster" (smaller arrival times) 
for small values of the inter- leg coupling t y . 

The clear difference in arrival times in the two regimes 
is striking. Some hints as to why it comes about can be 
found in the curve for t y /t x ~ 1. In the transition be- 
tween the two regimes (t y /t x ~ 1), a second peak, with 
a longer arrival time, becomes more prominent than the 
first peak. As t y increases, this second peak dominates 
and effectively marks the arrival of the doublon. In fact, 
this second peak arises as a result of the part of the dou- 
blon that travels through the lower leg instead of the 
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FIG. 6: (Color online) Excess double occupancy ANd(r) 
versus r for U — 10 and different values of t y . Inset: ANd at 
fixed times (r = and r = 6) versus t y . 

upper one for large t y , as shown in Figs. 14] and [5j 

A curious behavior is that, within botn regimes (t y < 
0.5 and t y > 2), the horizontal doublon speed increases 
with the coupling between the ladder legs. This is shown 
in the inset, where t~ 1 is plotted as a function of t y 
(circles). We notice that, within each regime, r" 1 (taken 
from the more prominent peak at each regime, except 
for t y — 1 where two values are used) does not depend 
linearly with t y , indicating that such increase does not 
arise from a simple rescaling of the "time unit" in terms 

of H/ty. 

In addition, the scaling of r" 1 with t y is different in 
each regime and therefore one might speculate if there is 
an extra energy scale dominating the behavior of t^ 1 at 
large t y . As it is known from previous studies P Hubbard 
two-leg ladders display a sizeable spin gap in the large t y 
regime. One possibility would be that the propagation 
of the doublon would be coupled to S = 1 spin excita- 
tions (where a singlet in a rung becomes a triplet), thus 
resulting in a scaling of the doublon speed with the spin 
gap. 

In order to check for such a connection between t^ 1 
and the spin gap A s , we have calculated A s for t y /t x > 2. 
The spin gap can be calculated with DMRG as 

A s =E(N/2 + l,N/2-l)-E(N/2,N/2) (3) 

where E(Nf,N±) is the energy of a state with N-f elec- 
trons with spin up and Ni electrons with spin down (to- 
tal N = Nf + N± electrons). We have performed these 
(static) calculations for a 2 x 32 (N = 64) ladder with 
U = 10 at half-filling. 

The results are shown in the inset of Fig. [7j Although 
the results do not show a clear scaling of the type "r^ 1 oc 
A s ", the dependence of r^ 1 with t y is markedly different 
in the strong spin-gap regime as compared to "small t y v 
one: the rate of increase of the doublon speed is lower for 
t y > 2t x , where a clear spin gap is present. This indicates 



FIG. 7: (Color online) Double occupancy at a site at the end 
of the ladder for (a) t y < 1 and (b) t y > 1 (no shift added). 
The position of the peak marks the "arrival time" r a of the 
doublon. Inset: inverse doublon arrival time r^ 1 (circles) and 
spin-gap calculated for a 2 x 32 ladder (squares) vs t v . Error 
bars are given by At/t 2 and indicate the uncertainty in the 
determination of the arrival time given the time step At. 

that the propagation of the doublon might be coupled to 
magnon-like spin excitations, although this connection is 
weak. This is also consistent with the rather weak decay 
shown in Fig. [6] even in the large t y regime. 

IV. CONCLUSIONS 

In this work, we have investigated the dynamics of 
holon-doublon pairs in Hubbard two-leg ladders using 
the time-dependent DMRG method. The geometry of 
the ladder brings some interesting qualitative changes as 
compared to the case of the chains. A telling example 
is a "non-local" character of the double occupancy: the 
creation of a doublon in a site at the upper leg reduces 
the double occupancy in the bottom site coupled to it. 
This, in fact, reduces the total double occupancy as the 
inter-leg coupling increases. 

In addition, we have found important qualitative dif- 
ferences in the dynamics of the excitations depending on 
the relative coupling between the ladder legs. For weak 
inter-leg coupling (t y /t x <C 1), results are qualit atively 
similar to the strictly one-dimensional caseP^M How- 
ever, the results show a strong downward shift in the 
doublon horizontal speed as one crosses over from weak 
(ty/tx <C 1) to strong (t y /t x >• 1) inter- leg coupling. We 
attribute this shift to different propagation "paths" for 
the doublon: for small t y , the doublon propagates mainly 
through the sites of a single leg. For large t y , the prop- 
agation involves also the vertical direction (through the 
ladder rungs rather than single sites), effectively "slowing 
down" the doublon. 

Although in both cases the "horizontal" speed of the 
doublon increases as the vertical coupling increases, the 
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scaling of the doublon speed with t y is different, being 
less pronounced in the strong inter-leg coupling regime. 
Our results indicate that the presence of a spin gap qual- 
itatively changes such scaling, indicating a possible (al- 
though weak) coupling between the doublon and spin ex- 
citations. 
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Appendix A: Accuracy during time evolution 

We have found that the accumulation of numerical er- 
rors during time evolution for the two-leg Hubbard ladder 
depends strongly on the ratio t y /t x . In addition to the 
usual criteria for detecting loss of accuracy during time 
evolution (namely, an exponential increase of the trunca- 
tion error for a given number m of kept sates, leading to 
a "runaway time"'^) , we have also monitored the energy 
expectation value E(t) = (*(t)|F|*(t))/<*(t)|*(t)) 
with increasing r as a way to probe the accuracy of the 
results for longer time scales. Since the Hamiltonian in 
Eq. [T] is time- independent, the expectation is that this 
quantity should be constant (i.e., E(t)/E(Q) — 1) and 
thus deviations from E(t)/E(0) — 1 can be used as a 
gauge of the accumulation of numerical errors during the 
time evolution of the system. 

The results for E(t) /E(0) are shown in Fig. [8] As time 
progresses, the accumulated errors during time evolution 



cause E(t)/E(0) to deviate from 1. This "degradation" 
in energy depends strongly on the inter-leg coupling t y : 
for t y in the ranges < t y < t x and t y > 2t x , results 
are within 5% accuracy up t ~ 6 using m = 400 states 
during time evolution and time step At = 0.2. For t y in 
the range 1.5 — 2, the degradation is stronger, although an 
accuracy of 10% (which we consider to be acceptable) up 
to t ~ 4 can be achieved. For t y > 2, the results improve 
again such that t y = 5 and t y = have comparable 
behavior. 

Keeping more states (say, m — 600) during time evo- 
lution brings only a slight improvement, as illustrated in 
the inset for t v /t x = 1 and we find that this increase does 
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FIG. 8: (Color online) Normalized average energy E(t)/ E(0) 
as a function of r for U — 10 with m — 400 kept states 
and different values of t y . Dashed and dotted horizontal lines 
mark the 5% and 10% "confidence ranges" respectively. Inset: 
E(t) for t y — 1 and different number of kept states during 
time evolution. 



not modify the trends and/or numerical values apprecia- 
bly: the results for the time-dependence of observables 
such as double occupancy are essentially unchanged up to 
7^6. For this reason, the calculations presented in this 
paper were performed keeping up to m = 400 states dur- 
ing time evolution, as increasing the number of retained 
states to m = 600 does not justify the considerable extra 
computational cost. 
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